home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / RANDCOMB.C < prev    next >
C/C++ Source or Header  |  1995-09-07  |  6KB  |  229 lines
  1. /*    rand_com[b].c
  2.     combination multiplicative random number generator
  3.         subtract two random numbers modulo moduli1-1 and shuffle
  4.     see
  5.         L'Ecuyer, Comm. of the ACM 1988 vol. 31
  6.         Numerical Recipes in C, 2nd edition, pp.  278-86
  7.         shuffling -- Knuth, vol. II
  8.     Copyright (c) 1994, 1995 by Gerald P. Dwyer, Jr.                */
  9.  
  10.  
  11. /* ------------------- */
  12. /* FUNCTION PROTOTYPES */
  13. /* ------------------- */
  14. # undef F
  15. # if defined(__STDC__) || defined(__PROTO__)
  16. #    define  F( P )  P
  17. # else
  18. #    define  F( P )  ()
  19. # endif
  20.  
  21. /* INDENT OFF */
  22. extern    long    genr_rand F((long, long, long, long, long));
  23. extern    long    genr_rand_diff F((void));
  24. extern    long    get_init_rand F((long));
  25. extern    void    init_rand_comb F((long, long));
  26. extern    long    rand_comb F((void));
  27.  
  28. # undef F
  29. /* INDENT ON */
  30.  
  31. #include    <time.h>
  32. #include    <stdlib.h>
  33. #include    <float.h>
  34. #include    <assert.h>
  35.  
  36. #define    TESTING
  37. # undef TESTING
  38.  
  39. #define    MOD1        2147483563L    /* first generator modulus */
  40. #define    MULT1        40014L        /* multiplier */
  41. /* ------------------------------------------- */
  42. /* modulus = multiplier * quotient + remainder */
  43. /* ------------------------------------------- */
  44. #define    Q1        53668L    /* quotient = [modulus/multiplier] */
  45. #define    R1        12211L    /* remainder */
  46.  
  47. /* second generator */
  48. #define    MOD2        2147483399L
  49. #define    MULT2        40692L
  50. #define    Q2        52774L
  51. #define    R2        3791L
  52.  
  53. #define    MOD_COMB    (MOD1-1)
  54.  
  55. #define    MIN_VALUE1    1
  56. #define    MAX_VALUE1    (MOD1-1)
  57. #define    MIN_VALUE2    1
  58. #define    MAX_VALUE2    (MOD2-1)
  59. #define    MAX_VALUE    ( (MOD1<MOD2) ? MAX_VALUE1 : MAX_VALUE2)
  60.  
  61. #define    GENR1(init_rand)    genr_rand(MULT1, init_rand, MOD1, Q1, R1)
  62. #define    GENR2(init_rand)    genr_rand(MULT2, init_rand, MOD2, Q2, R2)
  63.  
  64. #define    IMPOSSIBLE_RAND    (-1)
  65. #define    STARTUP_RANDS    16        /* throw away initial draws */
  66. #define    DIM_RAND    150        /* size of array shuffled */
  67.  
  68.  
  69. static long rand1, rand2 ;
  70. static long rand_num = IMPOSSIBLE_RAND ;
  71. static long rands[DIM_RAND] ;
  72.  
  73.  
  74. /* initialize generators with seeds
  75.     fill rands[] with initial values  */
  76. void init_rand_comb(long seed1, long seed2)
  77. {
  78.     extern long rand1, rand2 ;
  79.     extern long rand_num ;
  80.     extern long rands[] ;
  81.     int i ;
  82.  
  83.     if (seed1 <= 0 || seed1 > MAX_VALUE1)
  84.         seed1 = get_init_rand(MAX_VALUE1) ;
  85.     if (seed2 <= 0 || seed2 > MAX_VALUE2)
  86.         seed2 = get_init_rand(MAX_VALUE2) ;
  87.  
  88.                         /* seed the routine */
  89.     rand1 = seed1 ;
  90.     rand2 = seed2 ;
  91.  
  92.     genr_rand_diff() ;
  93.  
  94.     for (i = 1; i < STARTUP_RANDS; i++)    /* throw some away */
  95.         genr_rand_diff() ;
  96.                     /* fill the array of randum number values */
  97.     for (i = 0; i < DIM_RAND; i++)
  98.         rands[i] = genr_rand_diff() ;
  99.                     /* initialize rand_num for shuffling */
  100.     rand_num = rands[DIM_RAND-1] ;
  101. }
  102.  
  103.  
  104. /* get a long initial seed for generator
  105.     assumes that rand returns a short integer */
  106. long get_init_rand(long max_value)
  107. {
  108.     long seed ;
  109.  
  110.     srand((unsigned int)time(NULL)) ;    /* initialize system generator */
  111.     do {
  112.         seed  = ((long)rand())*rand() ;
  113.         seed += ((long)rand())*rand() ;
  114.     } while (seed > max_value) ;
  115.  
  116.     assert(seed > 0) ;
  117.  
  118.     return seed ;
  119. }
  120.  
  121.  
  122. /* generate the difference between random numbers
  123.     assumes    0 <  rand1    <  MOD1
  124.             0 <  rand2    <  MOD2
  125.     output       1 <= rand_num <= MOD_COMB */
  126. long genr_rand_diff(void)
  127. {
  128.     extern long rand1, rand2 ;
  129.     long rand_new ;
  130.  
  131.     rand1 = GENR1(rand1) ;
  132.     rand2 = GENR2(rand2) ;
  133.     rand_new = rand1 - rand2 ;
  134.     if (rand_new <= 0)
  135.         rand_new += MOD_COMB ;
  136.  
  137.     assert(rand_new >= 1 && rand_new <= MOD_COMB) ;
  138.  
  139.     return rand_new ;
  140. }
  141.  
  142.  
  143. /* generate the next value in sequence from generator
  144.     uses approximate factoring
  145.     residue = (a * x) mod modulus
  146.             = a*x - [(a*x)/modulus]*modulus
  147.     where
  148.         [(a*x)/modulus] = integer less than or equal to (a*x)/modulus
  149.     approximate factoring avoids overflow associated with a*x and
  150.         uses equivalence of above with
  151.     residue = a * (x - q * k) - r * k + (k-k1) * modulus
  152.     where
  153.         modulus = a * q + r
  154.         q  = [modulus/a]
  155.         k  = [x/q]        (= [ax/aq])
  156.         k1 = [a*x/m]
  157.     assumes
  158.         a, m > 0
  159.         0 < init_rand < modulus
  160.         a * a <= modulus
  161.         [a*x/a*q]-[a*x/modulus] <= 1
  162.             (for only one addition of modulus below) */
  163. long    genr_rand(long a, long x, long modulus, long q, long r)
  164. {
  165.     long k, residue ;
  166.  
  167.     k = x / q ;
  168.     residue = a * (x - q * k) - r * k ;
  169.     if (residue < 0)
  170.         residue += modulus ;
  171.  
  172.     assert(residue >= 1 && residue <= modulus-1) ;
  173.  
  174.     return residue ;
  175. }
  176.  
  177. /* get a random number from rands[] and replace it*/
  178. long rand_comb(void)
  179. {
  180.     extern long rand_num ;
  181.     extern long rands[] ;
  182.     int i ;
  183.                 /* if not initialized, do it now */
  184.     if (rand_num == IMPOSSIBLE_RAND) {
  185.         rand_num = 1 ;
  186.         init_rand_comb(rand_num, rand_num) ;
  187.     }
  188.  
  189.     /* rand_num from previous call determines next rand_num from rands[] */
  190.     i = (int) (((double)DIM_RAND*rand_num)/(double)(MAX_VALUE)) ;
  191.     rand_num = rands[i] ;
  192.  
  193.         /* replace random number used */
  194.     rands[i] = genr_rand_diff() ;
  195.  
  196.     return rand_num ;
  197. }
  198.  
  199. /* generates a value on (0,1) with mean of .5
  200.     range of values is [1/(MAX_VALUE+1), MAX_VALUE/(MAX_VALUE+1)]
  201.     to get [0,1], use (double)(rand_comb()-1)/(double)(MAX_VALUE-1) */
  202. double
  203. rand_rect_comb(void)
  204. {
  205.     return(double) rand_comb() / (double) (MAX_VALUE + 1);
  206. }
  207. #if    defined(TESTING)
  208. /* Test the generator */
  209. #include     <stdio.h>
  210. #define    EXP_VAL        804307721L
  211. void main(void)
  212. {
  213.     long seed1=1, seed2=1 ;
  214.     int i ;
  215.  
  216.     init_rand_comb(seed1, seed2) ;
  217.     printf("Seeds for random number generator are %ld   %ld\n",
  218.             seed1, seed2) ;
  219.     i = STARTUP_RANDS + DIM_RAND ;
  220.     do {
  221.         rand_comb() ;
  222.         i++ ;
  223.     } while (i < 9999) ;
  224.  
  225.     printf("On draw 10000, random number should be %ld\n", EXP_VAL) ;
  226.     printf("On draw %d, random number is        %ld\n", i+1, rand_comb()) ;
  227. }
  228. #endif /* TESTING */
  229.